Melting of a colloidal crystal 
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A melting transition for a system of hard spheres interacting by a repulsive Yukawa potential of 
DLVO form is studied. To find the location of the phase boundary, we propose a simple theory to 
calculate the free energies for the coexisting liquid and solid. The free energy for the liquid phase 
is approximated by a virial expansion. The free energy of the crystalline phase is calculated in the 
spirit of the Lenard-Jonnes and Devonshire (LJD) theory. The phase boundary is found by equating 
the pressures and chemical potentials of the coexisting phases. When the approximation leading to 
the equation of state for the liquid breakes down, the first order transition line is also obtained by 
applying the Lindemann criterion to the solid phase. Our results are then compared with the Monte 
Carlo simulations. 
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1. Introduction 

The asymmetric polyelectrolyte solutions provide a se- 
vere challenge, as well as, a true testing ground for a 
variety of statistical mechanics theories. Unlike the case 
of simple symmetric electrolytes [0, the phase structure 
of which is reasonably well understood, it is fair to say 
that our understanding of charge stabilized colloids is 
very far from complete. Even such basic questions as 
what is the form of interactions between strongly charged 
colloidal particles still remains controversial One 
thing, however, seems to be well accepted by both the ex- 
perimentalists and theorists, when the concentration of 
colloid is sufficiently large it freezes into a crystal [0-0 . 
The exact mechanism leading to the transition is still 
not fully understood [|l[ . The need to find a resolution to 
this dilemma is a particularly acute in view of industrial 
applications of charge stabilized colloids in the design of 
new materials, such as water soluble paints. 

The complete theory of charged colloids would have 
to include, in addition to polyions, the presence of neu- 
tralizing counterions. There have been a number of at- 
tempts to study this full system, however, due to se- 
vere mathematical complications, no attempt has been 
made so far to explore the liquid-solid phase transi- 
tion in highly asymmetric electrolytes [p^-[T7|. Instead 
most approaches rely on using an effective, counterion 
screened, interaction between the polyions of the DLVO 
form [^ 19|, while completely ignoring the contribution 
of counterions to the free energy. Although there have 
been some attempts to compare the phase diagrams de- 
rived on the basis of this "reductionist" approach directly 
with experiment, we should stress that there is no rea- 
son why this type of comparison should at all make sense 

At a first order phase transition the pressures and 
the chemical potentials of the coexisting phases must be 
equal. Since the biggest contribution to both pressure 
and chemical potential comes from the entropic motion 
of the counterions and the interactions between polyions 
and counterions are strongly coupled, any theory which 




does not explicitly take into account the presence of coun- 
terions is fatally flawed. Nevertheless, in the absence of 
a complete theory, the reductionist approach has an im- 
portant role as a testing ground of a variety of theories 
of the liquid-solid transition. However, the true merit 
of each theory should be measured by how it compares 
to the Monte Carlo (MC) simulation of the equivalent 
reduced model. 

Traditionally the theories of liquid-solid transition can 
be divided into two classes: those based on the density 
functional methods [^ and those based on the modifi- 
cations of the Lcnnard- Jones and Devonshire | ^ , p^ cell 
theory of liquid. Both methods have their strengths and 
weaknesses. 

In the case of simple electrolytes the cell model has 
proven to be in a closer agreement with MC [E3 
the early work of Fuoss, Katchalsky and Lifson 
cell theories formed the foundation on which most of the 
studies of polyelectrolyte solutions were based. The cell 
models are often indiscriminantly applied to both the 
solid and liquid phases. In this respect we would like to 
emphasize that the cell theories entail an underlying pe- 
riodicity which is valid only in the solid phase. From this 
perspective our understanding of a crystalline state of 
colloidal suspensions is quite satisfactory and agree 
with Tejero et al. approach for point particles [|l3|. What 
is lacking is an equally good treatment of a disordered 
state of suspension. Some preliminary approaches in this 
direction were recently reported in literature |p^ . In 
this paper we shall ignore many of the subtleties of the 
true colloidal suspensions, and confine our attention to 
a model system of hard spheres interacting by a Yukawa 
potential of the DLVO form. Furthermore we shall re- 
strict our attention to the strong screening limit, high 
salt concentration. In this case the repulsive interactions 
between the equally charged colloids will be sufficiently 
weak. Thus the virial expansion for, say, the pressure 
can be truncated at a leading order, allowing us to write 
an equation of state of a van der Waals form. Unfortu- 
nately, as the screening is reduced the truncation of the 
virial expansion will no longer be valid and an alterna- 
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tive method must be used. The crystalhne state of the 
suspension will be modeled by a modified Lennard-Jones 
and Devonshire theory. 

2. The model 

Our system consist of identical particles of charge —Zq 
and diameter a. The solvent is modelled as an uniform 
medium of a dielectric constant D. In the spirit of re- 
ductionist approach the presence of the counterions is 
neglected beyond their contribution to the screening of 
the interactions between the colloidal particles. We shall 
accept that this effective interaction is of a DLVO form 

n 

where k is an inverse Debye length, which sets the scale 
for the range of electrostatic interactions. It is convenient 
for later purposes to write the energy W in the form 



W = e (b(x) 



where 



X = r/a, e . 



{Zq? 



eff 



Da 



and the dimensionless potential (t>{x) is 
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In the limit Kcr ^ 1 , (j){x) is extremely short ranged. In 
this case all the thermodynamic properties of our system 
will be due to the hard-core repulsion between colloids. 
On the other hand, if kct <C 1 the thermodynamics will 
be dominated by the electrostatic interactions. 

Our aim will be to find the location of the liquid-solid 
coexistence region for the two limiting cases. To this end 
we will need the free energies for both liquid and solid 
phases. We shall now proceed to construct these free 
energies. 

3. The fluid phase (kct > 1) 

The Helmholtz free energy for the fluid phase is con- 
structed as a sum of terms, starting with an ideal gas 
contribution 

/3Ff = iV[log(pA3) - 1] , (5) 
where the thermal de Broglie wavelength is, 

A 



(27rmfcT)i/2 



The mass of each polyion is denoted by m, h is the Planck 
constant, and k is the Boltzmann constant. The hard- 
sphere repulsion between the polyions can be included 



through the Carnahan-Starling contribution to the free 
energy 
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where p* = pa^ is the reduced density. 

In the limit na ^ 1 the virial expansion for thermody- 
namic functions will be composed of rapidly decreasing 
terms. In this case it will be sufficient to approximate the 
electrostatic contribution to the free energy by a second 
virial term, 



N r 



Substituting (g) we find 

^ 2np*N '■+°" 



dx x^ (j){x) , 
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where the reduced temperature is T* = l//3e = 
Da/l3{Zq)1j:j:. The integral may be easily carried out, 
leading to 



(3Ff 
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The full free energy is Ff = F'f + Ff^ + Fj. It should 
be remembered that this expression is only valid for na ^ 
1; beyond this limit additional virial terms will contribute 

toF;. 

4. The solid phase 

We shall use the following mean field picture of a col- 
loidal crystal. Each polyion is allowed to move in a cage 
whose boundaries are defined by its nearest neighbors. 
However, the particles do not move freely since they in- 
teract electrostatically with the neighboring polyions. In 
the limit kct 3> 1 we shall consider the electrostatic in- 
teractions with the polyions of the first shell, although in 
principle other shells can be included reasonably straight- 
forwardly. 

The partition function can be written [ETl, 
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exp i-NpUo/2) 



(10) 



where Vf is the "free" or "effective" volume, available to 
the center of a polyion inside a cell. Note that no factors 
of A'^! are present, since in principle each particle of the 
lattice can be labeled. The term NUq/2 is the energy 
of the lattice when all the particles are located at their 
equilibrium positions. 

The free energy of the crystal is 



/3F, = ^/3C/„-A^ln(^) 



(11) 
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A final approximation of LJD theory is to simplify the 
geometry of the cell, by taking it to be a sphere with 
a surface charge Ge = ni{Zq)eff/4:TrR'^, where R = 
a{pcp/ py^^ is the distance between the nearest neigh- 
bors, Pep is the "close-packing" density for a particular 
crystalline form, and ni is the number of nearest neigh- 
bors. This is refered to as a "smearing procedure" used 
elsewhere (as in p8|). 

The free volume Vf is then simplified to 



An 



exp h/9(t/(r) - U{Q))ydr , (12) 



and U (r) denotes the energy of a particle inside the cell. 

In our case we shall assume that the solid state corre- 
sponds to an fee structure with p*^ = v^, since this is 
the transition found in the pure hard-sphere system . 

In order to find the mean field potential exerted on each 
polyion due to the neighboring polyions we must, in prin- 
ciple, integrate the effective potential (|^) over the surface 
of the sphere. This will require performing some suffi- 
ciently messy angular integrals (when the central polyion 
is displaced from the center of the cell). We can avoid 
this by observing that the DLVO potential satisfies the 
Helmholtz equation. 



(13) 



The required solution is unique, provided we consider 
the appropriate boundary conditions for the potential (j). 
These are, continuity at r = i?, finitencss inside the cell, 
vanishing value at infinity, and the discontinuity of the 
normal component of the electric field at r = i?, due to 
the presence of a surface charge. We find 
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The above integral is a function of the density and the 
temperature, and will be performed numerically. Given 
the free energy, the pressure, p = —dF/dV, and the 
chemical potential, p — dF/dN, of the liquid and solid 
phases can be calculated by a straightforward differen- 
tiation. At transition, the pressures and the chemical 
potentials of the coexisting phases are equal: 



PsiPs) = PfiPf) , f^siPs) = M/(P/) 



(19) 



Solving these equations numerically the boundary of the 
coexisting region can be obtained. The coexistence curve 
T* versus p* for three values of kct is ploted on Fig.l. 
The transition when the only interaction is due to the 
hard-spheres repulsion is represented by circles. We see 
that for large values of na our curve tends to these 
points, the electrostatic interaction being dominated by 
the hard sphere contribution. We also observe that in the 
limit of hi gh temperatures the electrostatic inter &tion 
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The free volume can now be obtained by performing 
the radial integral (^. We find 



Vf Aire" 



h , 



where 



sinhx 



exp — ai 



and ai is given by 



ai 



-kR 



kR 



(16) 



dx , X = Kr , (17) 



(18) 



FIG. 1. Fluid-solid coexistence for the repulsive Yukawa 
system. Three values of kct are shown: (a)K(T — 5, {b)Ka = 10, 
{c)Ka = 100. The circles represent the hard-sphere transition. 

The whole curve may be compared with that obtained 
in the MC simulation of reference ||^, and this is done 
in Fig. 2 for na = 5. The agreement is not very good, 
although the shape of the figure is nearly identical. Once 
again we would like to remind the reader that our the- 
ory remains valid only for sufficiently large values of na. 
When the screening weakens, the van der Waals approxi- 
mation will become less reliable, as more and more terms 
will need to be included in the virial expansion fsot]. 
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FIG. 2. Fluid-solid coexistence for the repulsive Yukawa 
system for kct = 5. The dotted lines are the MC results of [7]. 
The squares represent the result for point particles in [10]. 
We have ploted the results in the units of [7], PX = e~""^/r*, 
to allow comparison. 

In the absence of a more reliable expression for the 
free energy for the hquid state, we are, therefore, in need 
of a drastically different approach. One such method 
which relies purely on the information obtained from the 
solid state is the so called Lindemann criterion js^ for 
melting, which states that the crystal will melt when the 
amplitude of its vibrations become sufficiently large. In 
order to compare our results with the MC simulations for 
point particles |0 we set na = 0. 

The energy of a particle inside the cell is in this case 
given by 



f3U{r) = an 



sinh (nr) 



where 
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(20) 



(21) 



The semiempirical Lindemann criterion states that 
the solid will melt when the ratio of the mean-squared 
displacement ((r^)) to the square of the interparticle dis- 
tance (a^) is equal or greater than a constant of order 



10 . The quantity 



-1/3 



is a geometry indepen- 



dent term, hence applicable to the solid and to a disor- 
dered phase as the fluid. This gives an equation for the 
melting line, 



(22) 



with c a constant with typical values 0.025, 0.030, 0.0361, 
for Lennard- Jones types of potentials. 



The mean-squared displacement 
through the canonical average. 



/Q^r2exp(-/3C/(r))dT/ 
/o''exp(-/3C/(r))dF 



is obtained 



(23) 



The integrals must be done numerically, and we write 



(r 



_ h 



where 



47r 



kR 



Ii = —r dx X exp — oi 



47r 



kR 



dx X exp 



sinh; 



sinh a; 



-au 



(24) 



(25) 



(26) 



and X = KT. 

Therefore we may draw the melting line for several 
values of c. The resulting curves are the the solid lines 
shown in Fig. 3 for fee lattice. For typical values of c the 
curves are located somewhat below the simulation data 
of reference ||l^ and the theoretical results of Tejero et al. 
]l3[ . The agreement is better for higher values of c. This 
can be understood on the following grounds. We have 
considered only the first shell of neighbors, hence the 
potential inside a cell is found to be less repulsive than it 
is in reality. The net effect is that the amplitude of the 
vibrations of a molecule around its mean location in the 
cell is overestimated. This increases the value of the ratio 
(r^)/a^, which is precisely the Lindemann constant. 

This might be clearly seeing if we consider another ap- 
proximation at this level, which is to solve the integrals 
above analytically. For this purpose we extend the supe- 
rior limit of integration to infinity, and substitute the fac- 
tor svahx/x by the first two terms of its series, 1 -I- 2:^/6. 
Of course this will introduce an error, but we restrict 
ourselves to the limit of small Ka^, namely low salinity. 
In this case this approximation is acceptable, and the re- 
sult may be compared with that obtained from the exact 
numerical calculation. In this limit we find, 
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(27) 



From the above expression, the Lindemann criterion, 
Eq.(21), may be used to draw the melting line. The re- 
sults are ploted as the dashed lines in Fig. 3 for fee lattice. 
We see that for small values of kos the curves are close 
to the exact ones. Even though this simplification gives 
a worst agreement with the simulation, it has the ad- 
vantage of allowing us to write an analytical expression 
for the melting temperature as a function of the density. 
From the equation of melting, (p2|), and the expression 
for the relative mean squared displacement, Eq. (p7[) , we 
obtain the melting temperature. 



Tp(Kas) 



cni e 



367r kR 



(28) 



From the above equation we can see that the transition 
temperature depends on the coefficient cni. Therefore 
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FIG. 3. Melting lines for the repulsive Yukawa system as 
obtained from Lindemann criterion. We have point parti- 
cles and the fee symmetry. Solid lines are the result of exact 
calculation, and dashed lines are obtained from the gaussian 
approximation (see text). The full circles are the simulation 
coexistence points of [10]. The values of the Lindemann con- 
stant were taken as 0.025 in line (a), 0.0361 in line (6), and 
0.09 in line (c). 

The Lindemann criterion may also be used for parti- 
cles of finite diameter, and we do this for kct — 5, for 
comparison with the MC results of |Q. In this case we 
have 



{r ) = -r , 

12 



where 



_ in 
J 1 = ^ 



Jo 



kR—kct 



dx X exp — ai 



sinha; 



X 



(29) 



(30) 



of the constant c, and the results are the dashed lines in 
Fig. 4. With the above expression for K^(r^), the equation 



for the melting, ( |22| ) , we obtain the melting temperature 
as for point particles, T* = T*Ka. Note the difference 
in the definitions of reduced temperatures for finite size 
and for poliit particles.' ' ' ' /' 

We note that once again the-^e,S'a^9'ement wit iqcom- 
puter experiment is found for c 4: 0.09./ Furthefmore the 
approximated form for the melti^gi'teniijaeralute *pr jduces 
a better fii- to MC, sugesting tlp.il the/cell model 
timates th; effects of confinemdnt. ' " 




FIG. 4. Melting lines for the repulsive Yukawa system as 
obtained from Lindemann criterion. The particles have finite 
size(K(j = 5) and the symmetry is fee. Solid lines are the 
result of exact calculation, and dashed lines are obtained from 
the gaussian approximation (see text). We have adopted the 
units of [7] to allow comparison. The dotted lines are the 
MC results of [7]. The squares represent the result for point 
particles in [10] . The values of the Lindemann constant were 
taken as 0.025 in line (a), 0.0361 in line (6), and 0.09 in line 
(c). 



5. Conclusions 



47r f^^-'^'^ ( sinha;\ , , 

12 = — J ax X exp y—ai j , (31) 

and Oi is defined in Eq.(p^). 

With the above results the melting line is constructed, 
and the results are shown as the solid lines in Fig. 4. The 
dotted line is the curve of [^ . Wc now consider a simpli- 
fied version as was done in the point particle case, namely 
to evaluate the integrals in the expression for (r'^) from 
to oo, and to replace smhx/x by 1 + x'^/6. The same 
range of validity is understood. We obtain as before 




but remember that ai is now given by Eg. (p^ . We may 
draw the corresponding melting lines for several values 



We have studied a melting transition in a system com- 
posed of hard spheres interacting by a repulsive Yukawa 
potential. A simple theory to calculate the free energies 
of the coexisting liquid and solid phases was proposed. 
Since our equation of state for the liquid is based on the 
virial expansion, it is expected to work well only in the 
limit of strong screening, since in this limit it is sufficient 
to keep the first virial term. However, if a more reliable 
equation for the liquid state can be found it should, in 
principle, be possible to extend our results over the full 
range of salinity. For the solid phase, we introduced a 
cell model based on the Lennard-Jones and Devonshire 
theory. From the free energy we derived expressions for 
the pressure and chemical potential. 

The coexistence curve was obtained by equating the 
pressures and the chemical potentials of both phases. We 
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compared our results with simulations in Fig. 2. Even- 
though the results are qualitatively similar, the actual 
value of the densities for the coexisting phases are quite 
different. This is most likely due to the breakedown in 
our approximation leading to the equation of state for 
the liquid. 

In the absence of a better liquid state theory, we have 
relied on the Lindemann criterion to estimate the loca- 
tion of the melting line for low salinity. With a suitable 
choice for the Lindemann parameter, a reasonable fit to 
the MC data was found. Unfortunately this parameter 
is significantly larger than is usual for the systems with 
short range interactions. 
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